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Abstract 

Hit-and-Run is known to be one of the best random sampling algorithms, its 
mixing time is polynomial in dimension. Nevertheless, in practice the number 
of steps required to achieve uniformly distributed samples is rather high. We 
propose new random walk algorithm based on billiard trajectories. Numerical 
experiments demonstrate much faster convergence to uniform distribution. 
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1. Introduction 

Generating points uniformly distributed in an arbitrary bounded region 
Q C W 1 (sampling) finds applications in many computational problems [1, 2]. 

Straightforward sampling techniques are usually based on one of three ap- 
proaches: rejection, transformation and composition. In rejection approach 
the region of interest Q is enclosed within the region with available uniform 
sampler B (usually a box or a ball). At the next step non belonging to Q 
samples are rejected. Suppose Q is an unit ball while bounded region B 
is a box [—1, l] n . Then for n = 2k we obtain the ratio of volumes of the 

box and the ball q = — — — - = — r , thus q =~ for n = 20 and we 

Vol(S) k\2 k 

should generate ~ 10 8 samples to have just a few in Q. For polytopes this 
ratio can be much smaller. The other way to exploit pseudo-random number 
generator for simple region B is to map B onto Q via smooth deterministic 
function with constant Jacobian. For instance, to obtain uniform samples in 
Q = {x : x T Ax < 1}, A being positive definite matrix, it suffices to take y 
uniform in the unit ball ||y||2 < 1 and transform them A- 1/2 y. Un- 

fortunately, such a transformation exists just for a limited class of regions. 
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In composition approach we partition Q for finite number of sets that can 
be efficiently sampled. Apart from narrow class of regions with available 
partition, Q is partitioned into finite union of simplices and the number of 
simplices makes the procedure computationally hard. 

Other sampling procedures use modern versions of Monte Carlo tech- 
nique, based on Markov Chain Monte Carlo (MCMC) approach jE 0|. For 
instance, recent efficient algorithms for volume computation based on random 
walks can be found in IE [g] . One of the most famous and effective algorithms 
of MCMC type is Hit-and-Run (HR) that was originally proposed by Turchin 
Q and independently by Smith [g]. The examples of Hit-and-Run applied to 
various control and optimization problems are provided in (EEE 11 1. Unfor- 
tunately, even for simple "bad" sets, such as level sets of ill-posed functions, 
HR techniques fail or at least are computationally inefficient. A variety of 
applications and drawbacks of existing techniques propose much room for 
improvement new sampling algorithms. For instance, there were attempts to 
exploit the approach, developed for interior-point methods of convex opti- 
mization [lij, and to combine it with MCMC algorithms. As a result Barrier 
Monte Carlo method fl3| generates random points that are preferable in 
comparison with standard Hit-and-Run. But the complexity of each itera- 
tion in general is high enough (the calculation of (V 2 -F(x)) , where F(x) 
is a barrier function of the set, is needed). Moreover such approach can not 
accelerate convergence for sets like simplices. 

In this paper we propose the new random walk algorithm, which is moti- 
vated by physical phenomena of a gas diffusing in a vessel. A particle of gas 
moves with constant speed until it meets a boundary of the vessel, then it re- 
flects (the angle of incidence equals the angle of reflection) and so on. When 
the particle hits another one, its direction and speed changes. In our simpli- 
fied model we assume that direction changes randomly while speed remains 
the same. Thus our model combines ideas of Hit-and-Run technique with 
use of billiard trajectories. There exist a vast literature on mathematical bil- 
liards, and many useful facts can be extracted from there 

Q EE HE QUI- 

Traditional theory addresses the behavior of one particular billiard trajec- 
tory in different billiard tables, their ergodic properties and the conditions 
for existence of periodic orbits. We extend billiard trajectories with random 
change of directions, this introduction of randomness enriches their ergodic 
properties. 

The paper is organized as follows. In Section 2 we present novel sampling 
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algorithm and prove that it produces asymptotically uniformly distributed 
samples in Q. Section 3 is dedicated to discussion of some properties of the 
Billiard Walk, implementation issues are discussed as well. Simulation of 
BW for particular test domains is presented in Section 4. Much attention is 
devoted to ability of BW to get out of the corner in comparison with HR. 
Here we consider just the most demonstrative types of geometry. In Section 
5 we briefly discuss possible applications of the algorithm. 

2. Algorithm 

Suppose there is a bounded closed connected set Q C K" and a point 
a; G Q. Our aim is to generate asymptotically uniform samples x l G Q, 
' 1 v. 

The brief description of Hit-and-Run algorithm is as follows. At every 
step HR generates a random direction uniformly over the unit sphere and 
chooses next point uniformly from the segment of the line in given direction 
in Q. 

New algorithm Billiard Walk (BW) generates a random direction uni- 
formly as Hit-and-Run. But the next point is chosen as the end of the 
billiard trajectory of length t. This length is chosen randomly: we assume 
that probability of collision with another particle is proportional to 5t for 
small time instances 5t, this validates the formula for I in algorithm below. 
The scheme of the method is given in Fig. [T] while the precise routine is as 
follows. 

Algorithm of Billiard Walk. 

1. Starting point x° G Int Q is given; i = 0, x = x . 

2. Generate the length of the trajectory £ = — rlog£, £ being uniform 
random in [0, 1], r is a specified parameter of the algorithm. 

3. Pick random direction d G M n uniformly distributed on the unit sphere 
(i.e., ef = C/IICH, £ = randn (n, 1) - the n- dimensional vector with nor- 
mally distributed components). Construct billiard trajectory starting 
at x l and with initial direction d = d l . When the trajectory meets a 
boundary with internal normal s, ||s|| = 1, the direction is changed as 

d d — 2(d, s)s. 

4. Calculate the end of the trajectory of length £. If the point with nons- 
mooth boundary is met or the number of reflections exceeds 10n go to 
step 3. 
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5. i — i + 1, take the end point as x t+1 and go to step 2. 




-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 



Figure 1: Billiard Walk. 

Theorem 1. Suppose Q is connected, bounded and open (or a closure of 
such set) set, the boundary of Q is piecewise smooth. Then the distribution 
of points x l sampled by BW algorithm tends to uniform on Q. 

Proof. First, the algorithm is well defined: with probability one x t+1 will be 
found for arbitrary x 1 G Int Q (i.e. billiard trajectory with arbitrary initial 
conditions can be extended for arbitrary finite length). Note that we restrict 
the number of reflections to avoid situations when the trajectory comes (in 
limit) to a point with no normal with nonzero probability and for fixed I the 
length of billiard trajectory remains < £ after any number of reflections. 

All the restrictions for Q are important. Connectedness guarantees that 
stating from any point we can reach any other point of Q. Boundedness is 
necessary to define uniform distribution on Q and to avoid trajectories going 
to infinity. Openness allows us to connect any two points with a tube of 
nonzero measure. Thus, there exists piecewise linear trajectory connecting 
arbitrary points. 



4 



In view of Theorem 1 in |8j it suffices to prove that p(y\x) > for all x, y £ 
IntQ and p(y\x) = p(x\y), where p(y\x) is probability density of transition 
from x to y. Inequality p(y\x) > is guaranteed because all the directions 
are possible, Q is connected and open, and probability of any length i is 
positive. Equality p(y\x) = p(x\y) (reversibility) holds due to reflection law: 
the angle of incidence equals the angle of reflection. Thus both conditions are 
satisfied and the distribution of points x l sampled by BW algorithm tends 
to uniform on Q. □ 

3. Discussion 

We start discussion from some implementation issues. 

3.1. Choice of t . 

To run the algorithm parameter r need to be specified. The value of r 
strongly influences the behavior of the method. When r is small enough BW 
becomes slower that HR, it behaves as a ball walk with radius r. Empirical 
observations show that fast convergence to uniform distribution is achieved 
for r ~ diamQ. 

3.2. Preliminary transformation of Q. 

If Q is "ill-shaped", sometimes it can be improved with its linear transfor- 
mation. For instance, if Q is a box Q = {x : \xi\ < <2j} with a; having large 
data scattering simple scaling transforms Q into a cube. Similar scaling can 
be done to transform an ellipsoid into a ball. In general case the following 
scaling can be helpful. Suppose Q has a barrier F(x) [l2j|; for instance, for the 
polytope Q = {x : (ai,x) < hi] this barrier is F{x) = — ^^log(6j — (a^x)). 

i 

Then it is an easy task to find x* — an approximate minimum of F(x). Dikin 
ellipsoid E = {x : (H(x — x*), (x — x*)) < 1}, H — V 2 F(x*) lies in Q and 
is a good approximation of the polytope Q. Hence we can calculate linear 
mapping T = if -1 / 2 ; generating directions d! = Td, where d is uniformly 
distributed on the unit sphere, we can strongly accelerate the convergence. 
However sometimes none of transformations can improve the shape of the 
set, a simplex is known to be the worst-case example. 
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3.3. Boundary oracle and normals 

Both HR and BW algorithms require the calculation of intersection of a 
straight line (defined by point x k and direction d of the trajectory) with the 
set Q. We call this segment Boundary Oracle (BO) that is the segment [t,t\, 
where 

t = maxit : x k + td E dQ\, t = minit : x k + td G dQ} 

t<o t>o 

(we suppose that Q is convex, otherwise the points of first intersection of 
straight line and boundary of Q are taken). In most applications finding BO 
is not a problem. For instance, if Q is a polytope 

Q = {x G R n : (a 1 , x) <bi, i = 1, . . . , m} 

then 

hi — (a 1 , x k ) 



(a 1 , d) 



i = l, . . . , m, 



maxi, t = mint. 
t<o t>o 



Numerous examples of BO for other sets QJfor instance, defined by Linear 



Matrix Inequalities) can be found in (9), [Io|, 11 

BW walks requires also calculation of normals s in boundary points. In 
most applications it is not hard, for instance, for a polytope s = ai, where % 
is the index, for which maximum or minimum in above formulas is achieved. 



3.4- Comparison of HR and BW 

Our goal in test examples below is to compare HR and BW. We use 
several tools for this purpose. Sometimes theoretical considerations can help 
to compare the number of iterations to quit from the corner. It is well 
known, that HR can require too many iterations to get out of the corner, see 
for instance estimates in jl9|. We will show that estimates for BW are much 
more optimistic for many particular examples. On the other hand, we use 
simulation for comparison as well. We exploit different tools to demonstrate 
that one sampling is closer to uniform than another. Sometimes graphical 
figures in 2D plane are quite evident. In other cases we demonstrate strong 
serial correlation in samples (for instance, distance to initial point increases 
too slowly). Finally, we use parametric partition of Q and compare the 
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number of empirical points compared with theoretical number for uniform 
distribution. 

To make final conclusions on comparison of both methods we should have 
in mind the following. Of course, computationally BW is harder than HR. It 
requires more BO calculations, each reflection at the boundary also requires 
extra calculations, and it is not obvious how to estimate the complexity of 
each iteration of BW. Nevertheless acceleration of convergence to uniform 
distribution often makes BW preferable if compared with HR. 



4. Test sets and simulation 

Some sets below are unbounded; they are considered just to analyze the 
behavior at a corner. We say that a trajectory quits the corner if it goes to 
infinity. 

4-1. Angle at a plane 

Let Q C R 2 be an angle equal a < ir. Then billiard trajectory inde- 
pendently of initial point and initial direction quits Q after no more than 
N* = [7i/<y\ reflections, here \_a\ stands for smallest integer > a. It can 
be proved this way: if we reflect our angle N times around its side, billiard 
trajectory will be a straight line. It can not intersect any straight line twice. 

For HR we quit Q with probability 1 — (1 —a/ir) N after N iterations. For 
N = N* large we quit Q with probability 1 — 1/ e = 0.63 after iV* iterations, 
while for BW we do it w.p.l. 

It is of interest to estimate an average number of reflections (over random 

initial directions). Consider the triangle Q = {x E M. 2 : \xi\ < atan— ,X2 < 

1} with angle a. Let BW trajectories start at [0;0.1] and calculate the 
number of reflections until the trajectory reaches the line xi = \. For a = 7r/4 
25 trajectories are plotted in Fig. [2j 

The results for 5000 runs and various a are given in Table 1. Note that 
average number of reflections equals N* /2. 



4-2. Multidimensional case - polyhedral cone Q 

For polyhedral cone there exists number M independent of initial data 
such that billiard trajectory quits Q after no more than M reflections (see 
17J, also 14| . Theorem 7.17). However M depends on geometry of Q. If 
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Figure 2: 25 trajectories reaching line X2 — 1 starting from [0;0.1]. 



a 


Mean 


Std 


tt/2 


2.28 


0.87 


vr/4 


3.08 


1.3 


7I-/10 


5.94 


2.93 


tt/50 


25.08 


14.46 



Table 1: Mean and standard deviation for the number of reflections required to quit the 
corner with angle a. 



M is large (M > lOn) then BW algorithm sometimes returns to the initial 
point. However it is possible to prove, that BW is well defined w.p.l. 

4.3. Orthant Q = {x G R n : x > 0} 

It is easy to note that billiard trajectory independently of initial point 
and initial direction quits Q after no more than n reflections. Indeed, if d is 
direction of trajectory, I = {i : di < 0} then at each reflection / decreases, 
and after < n reflections 7 = 0. 

HR trajectory quits Q with probability 1 — (1 — 2 n_1 ) after a single itera- 
tion, thus it requires approximately 2 n_1 iterations to quit Q with probability 
1 — 1/e = 0.63. Hence BW is much more effective than HR for this case. 
Simulations for a cube (see below) confirm this statement. 
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4-4- Convex corner 

All these results show that a polyhedral corner is not a problem for BW in 
contrast with HR, where distance of the initial point to the boundary plays 
a significant role. The results can be extended for curvilinear corners with 
nondegenerate linear approximation, i.e. if linear approximation of a corner 
is a polyhedral cone with nonempty interior. 



4-5. Concave corner 



It is known 17] that concave corners can be attraction points for billiard 



trajectories. Consider typical set 

Q = {x e R 2 : H^Hoo < 1, \\x- ciiH > 1} (1) 

cij are vertices of HxH^ < 1. Then Q has 4 concave corners (cusps). Simula- 
tion of the single trajectory starting in [—0.001; 0] in direction [0; —1] is shown 
in Fig. [3l It requires 121 reflection to get out of the corner. For x = [— 10~ 4 ; 0] 
and the same direction it requires 670 reflections, for x = [— 10 -5 ; 0] the num- 
ber of reflections is 3802. But, in general, these "bad" directions are rare. Fig. 
H] (left) depicts 200 points for the set ([I]), the average number of reflections 
per point is 6. 




-1 -0.5 0.5 1 ' -1.S -1 -0.5 1 0.5 t 1.5 



Figure 3: Left: the trajectory for the set (TTJ) starting in [— 0.001; 0] in direction [0;— 1]. 
Right: zoom of the lower part. 



Consider slightly different concave angle 

Q = {x G R 2 : -x\ < x 2 < xf, Xi > 1} (2) 

with curvature tending to zero. It may happen that billiard trajectory comes 
(in limit) to a point with no normal with nonzero probability and for fixed 
I the length of billiard trajectory remains remains < I after any number of 
reflections. Indeed, start a trajectory at point x° = [0.9; e], e being small 
enough, fix i = 1, d — [—1; 0] and calculate the number of reflections needed 
to realize the trajectory. The results are shown on Table 2. 



e 


Number of reflections 


le-3 


746 


5e-4 


1851 


4e-4 


2480 


3e-4 


3617 


2e-4 


6158 


l.le-4 


13496 


1.01e-4 


>5e+6 



Table 2: Number of reflection required to realize the trajectory of length 1 for domain (J21 
starting from x° = [0.9; e] in direction d = [— 1;0]. 

As one can notice the number of reflections increases dramatically as the 
first coordinate of x° tends to zero and even for x\ = 10~ the trajectory 
becomes unreliazable. So to be on the safe side of situations like this we 
restrict the number of reflections in the algorithm. In Fig. H] (right) 200 
samples for domain fl2]) are depicted, average number of reflections is 5.2. 

4 ■ 6. Strip 

For a set like Q = {x E R 2 : < x 2 < I, \xi\ < M}, M being large 
enough HR and BW have different abilities to walk along X\. If we count 
average step per 1 boundary oracle application, BW is approximately 4 times 
faster. Indeed, HR takes 2 BO for 1 step and takes approximately the middle 
of a segment, while BW for 1 BO takes the segment height and its shift along 
x\ appears to be twice larger than for HR. 
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Figure 4: Left: 200 point for domain {TJ. Right: 200 point for domain ([2]). 
4.7. Cube 

For the unit cube Q = {x G IR n : < x < 1} (inequality understood 
component-wise) we can derive the next point of the BW algorithm explicitly. 

At the current point x for given I and d calculate ki = \_Xi + ld%\ ( \_x\ is 
maximal integer less than or equal to x) and walk to y: 

J Xi + Edi — ki, ki is even . _ 

Vi= \ l-( Xi + idi -ki), h is odd ' 1 = ' " " " ' n ' 

Of course there is no need to apply MCMC algorithms for random sampling 
in a cube: a code rand(n, 1) in MatLab generates uniformly distributed 
points. Moreover, the shape of a cube is so nice that distribution of HR 
points converges to uniform fast enough. Nevertheless serial correlation for 
these points is evident. In Fig. [5] we compare = E\ \x k — x°\ for n = 50 
averaged over 500 runs for two initial points x° = [1/2, . . . , 1/2] T (left) and 
x° = [1/n, . . . , ^/n\ T (right). Implementing BW we take r = \fn. One can 
see that serial correlation is much stronger for HR samples (black) than for 
BW samples (blue), and even 100 iterations are not enough for HR algorithm 
to quit the corner. 
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Figure 5: n = 50. £711 X X 1 1 qq cLVGI"cl ged over 500 runs for a; = [1/2, ... , 1/2] T (left) and 
x° = [1/n, . . . , l/n} T (right) for HR samples (black) and BW samples (blue). 



4-8. Simplex 

The next test set is standard n-dimensional simplex 

Q = {xj > 0, },Xj = l,i = 0, 1, . . . ,n}. 

Simplex is a set, containing many corners, while the geometry of simplex 
can't be improved by any affine transformation. We know that for HR walk 
it takes a lot of iterations to get out of a corner, thus it is interesting to 
compare HR and BW. 

Smooth boundary of Q is specified by points B = {x G 
for one k} and internal normal of unit length at such points is 

T 



pn+1 



: x k 



1 



n(n + 1) 



■1,... 



n 



Edge length is for every dimension n and we take r = \/2. Fig. [6] shows 
300 points generated by HR (black) and BW (blue) for standard 2-simplex. 

We see that for n = 2 samples look uniformly distributed for both algo- 
rithms. To judge about the uniformity more rigorously in multidimensional 
case, consider the sequence of enclosed simplices S a = {x e M n+1 : X{ > 



a 



E 



.1 i 



1}, < a < 



n 



For a = we have the initial simplex, for 
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Figure 6: 300 points generated by HR (black) and BW (blue) points for standard 2-simplex. 



a = simplex S a contains one point. Let f(a) be a portion of points 

71 + 1 

contained in S a , and denote /(a) = vol S a / vol Sq = (1 — (n + l)a) n . Fig. [7] 
shows /(or) for n = 50, JV = 300,2° = {l/(n + 1), . . . l/(n + 1)}. Red line 
corresponds to uniformly distributed points, black line describes the distri- 
bution for HR points and blue line for BW points. We conclude that for BW 
samples empirical values of f(a) are much closer to mean value f(a) than 
for HR samples. 

Another advantage of BW — its ability to quit a corner fast enough. We 
take a specific starting point close to the corner x° ' ' 
0.1. 

Fig. [8] depicts the distance between x 1 and the vertex of the simplex 
c = [1,0,..., 0] T . HR is unable to get out of the corner while the behavior of 
BW points after a few iterations looks like for uniformly distributed points. 



l-e 

n nJ 
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Figure 7: Portion of points contained in S a for uniformly distributed points (red), HR 
(black) and BW (blue), n = 50, 300 points. 
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Figure 8: Distance from the corner [1,0,..., 0] T for the first 100 points of walks in 50- 



simplex for different starting points. x° = 1 — i 
tributed points (red), HR (black) and BW (blue). 
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4-9. Toroid 

Both algorithms HR and BW are applicable for nonconvex sets. Consider 
the toroid formed by a n-dimensional ball of radius r with its center rotating 
in a circle in (xi, X2)-plane: 

Q = {x E R n : \\x-c x \\ < r}, (3) 

where c xi = X% , % = 1, 2, = 0, z > 2. 
Vxf + xi. 

Fig. [9] demonstrates = 1000 samples (projected onto (x\, a^-plane) 
for the set ([3]) with r = 1/3 of dimension 10 (x\, X2)-plane. HR points are 
plotted with black dots, BW points with blue. 




Figure 9: (xi, X2)-projection for HR points (black) and BW points (blue) for toroid |3|). 
n = 10, N = 10 3 . 

It can be easily seen that angles of BW points are much more uniformly 
distributed than for HR points, which remain in the neighborhood of the 
initial point. 

5. Applications 

In this paper we do not address numerous applications of new version of 
random sampling. We can mention just few of them — global optimization 
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(in particular, concave programming), control problems, robustness issues, 
numerical integration, calculation of volume and of center of gravity and so 
on, see, for instance, our previous papers [Io|, 11, 13 . We plan to consider 
these applications in future works. 
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